Functions from these packages will be used throughout this document:
[R code]
library(conflicted) # check for conflicting function definitions# library(printr) # inserts help-file output into markdown outputlibrary(rmarkdown) # Convert R Markdown documents into a variety of formats.library(pander) # format tables for markdownlibrary(ggplot2) # graphicslibrary(ggfortify) # help with graphicslibrary(dplyr) # manipulate datalibrary(tibble) # `tibble`s extend `data.frame`slibrary(magrittr) # `%>%` and other additional piping toolslibrary(haven) # import Stata fileslibrary(knitr) # format R output for markdownlibrary(tidyr) # Tools to help to create tidy datalibrary(plotly) # interactive graphicslibrary(dobson) # datasets from Dobson and Barnett 2018library(parameters) # format model output tables for markdownlibrary(haven) # import Stata fileslibrary(latex2exp) # use LaTeX in R code (for figures and tables)library(fs) # filesystem path manipulationslibrary(survival) # survival analysislibrary(survminer) # survival analysis graphicslibrary(KMsurv) # datasets from Klein and Moeschbergerlibrary(parameters) # format model output tables forlibrary(webshot2) # convert interactive content to static for pdflibrary(forcats) # functions for categorical variables ("factors")library(stringr) # functions for dealing with stringslibrary(lubridate) # functions for dealing with dates and times
Here are some R settings I use in this document:
[R code]
rm(list =ls()) # delete any data that's already loaded into Rconflicts_prefer(dplyr::filter)ggplot2::theme_set( ggplot2::theme_bw() +# ggplot2::labs(col = "") + ggplot2::theme(legend.position ="bottom",text = ggplot2::element_text(size =12, family ="serif")))knitr::opts_chunk$set(message =FALSE)options('digits'=6)panderOptions("big.mark", ",")pander::panderOptions("table.emphasize.rownames", FALSE)pander::panderOptions("table.split.table", Inf)conflicts_prefer(dplyr::filter) # use the `filter()` function from dplyr() by defaultlegend_text_size =9run_graphs =TRUE
1.1 Overview
Time-to-event outcomes
Survival analysis is a framework for modeling time-to-event outcomes. It is used in:
clinical trials, where the event is often death or recurrence of disease.
engineering reliability analysis, where the event is failure of a device or system.
insurance, particularly life insurance, where the event is death.
Note
The term survival analysis is a bit misleading. Survival outcomes can sometimes be analyzed using binomial models (logistic regression). Time-to-event models or survival time analysis might be a better name.
1.2 Time-to-event outcome distributions
Distributions of Time-to-Event Data
The distribution of event times is asymmetric and can be long-tailed, and starts at 0 (that is, \(P(T<0) = 0\)).
The base distribution is not normal, but exponential.
There are usually censored observations, which are ones in which the failure time is not observed.
Often, these are right-censored, meaning that we know that the event occurred after some known time \(t\), but we don’t know the actual event time, as when a patient is still alive at the end of the study.
Observations can also be left-censored, meaning we know the event has already happened at time \(t\), or interval-censored, meaning that we only know that the event happened between times \(t_1\) and \(t_2\).
Analysis is difficult if censoring is associated with treatment.
Right Censoring
Patients are in a clinical trial for cancer, some on a new treatment and some on standard of care.
Some patients in each group have died by the end of the study. We know the survival time (measured for example from time of diagnosis—each person on their own clock).
Patients still alive at the end of the study are right censored.
Patients who are lost to follow-up or withdraw from the study may be right-censored.
Left and Interval Censoring
An individual tests positive for HIV.
If the event is infection with HIV, then we only know that it has occurred before the testing time \(t\), so this is left censored.
If an individual has a negative HIV test at time \(t_1\) and a positive HIV test at time \(t_2\), then the infection event is interval censored.
1.3 Distribution functions for time-to-event variables
The Probability Density Function (PDF)
For a time-to-event variable \(T\) with a continuous distribution, the probability density function is defined as usual (see probability density function).
Example 1 (exponential distribution) Recall from Epi 202: the pdf of the exponential distribution family of models is:
Theorem 2 If \(A_t\) represents survival status at time \(t\), with \(A_t = 1\) denoting alive at time \(t\) and \(A_t = 0\) denoting deceased at time \(t\), then:
\[\text{S}(t) = \Pr(A_t=1) = \text{E}[A_t]\]
Theorem 3 If \(T\) is a nonnegative random variable, then:
Definition 3 (Incidence rate) Given a population of \(N\) individuals indexed by \(i\), each with their own hazard rate \({\lambda}_i(t)\), the incidence rate for that population is the mean hazard rate:
Theorem 4 (Incidence rate in a homogenous population) If a population of individuals indexed by \(i\) all have identical hazard rates \({\lambda}_i(t) = {\lambda}(t)\), then the incidence rate for that population is equal to the hazard rate:
\[\bar{{\lambda}}(t) = {\lambda}(t)\]
Lemma 1 (Joint probability of a variable with itself)\[p(T=t, T\ge t) = p(T=t)\]
Proof. Recall from Epi 202: if \(A\) and \(B\) are statistical events and \(A\subseteq B\), then \(p(A, B) = p(A)\). In particular, \(\{T=t\} \subseteq \{T\geq t\}\), so \(p(T=t, T\ge t) = p(T=t)\).
Theorem 5 (Hazard equals density over survival)\[{\lambda}(t)=\frac{f(t)}{\text{S}(t)}\]
As we will see below, \({\Lambda}(t)\) is tractable to estimate, and we can then derive an estimate of the hazard function using an approximate derivative of the estimated cumulative hazard.
Example 5 The cumulative hazard function for the exponential distribution with rate parameter \(\lambda\) is:
\[
{\Lambda}(t) = \mathbb{1}_{t \ge 0}\cdot \lambda t
\]
Figure 3 shows some examples of exponential cumulative hazard functions.
Figure 3: Examples of exponential cumulative hazard functions
Some Key Mathematical Relationships among Survival Concepts
If instead the event time \(T\) is censored and only known to be after time \(y\), then the likelihood of that censored observation is instead the survival function evaluated at the censoring time:
Exercise 3 (Construct a prediction interval) Suppose a cancer patient is predicted to have an expected (mean) lifetime of 7 years after diagnosis, and suppose the distribution is exponential.
Construct a 95% prediction interval for survival.
Tip
Use the quantiles of the exponential distribution.
Solution 1. If the mean is 7 years until death, then the rate parameter \(\lambda = 1/7\) events (deaths) per year.
The 0.025 quantile of the exponential distribution with \(\lambda = 1/7\) is qexp(p 0.025, rate = 1/7) = 0.177225 and the 0.975 quantile is qexp(p 0.975, rate = 1/7) = 25.822156, so the prediction interval is qexp(p c(.025, 0.975), rate = 1/7) = (0.177225, 25.822156).
Exercise 4 Graph the prediction interval as a function of the mean, for Gaussian (\(\sigma = 1\)), Binomial, Poisson, and Exponential.
Solution 2. Left to the reader for now.
Exercise 5 (Explain the results) Why do time-to-event models have such wide predictive intervals?
Tip
Consider the relationship between the mean, variance, and standard deviation of the exponential distribution, and contrast that relationship with the Poisson distribution and the Bernoulli distribution.
Solution 3. In the exponential distribution, variance is the square of the mean (hence SD is equal to mean); as opposed to Poisson, where variance was equal to the mean (and SD is the square-root of the mean), or the Bernoulli, where the variance is the mean minus the square of the mean (so the SD is smaller than the square-root of the mean).
Estimating \(\lambda\)
Suppose we have \(m\) exponential survival times of \(t_1, t_2,\ldots,t_m\) and \(k\) right-censored values at \(u_1,u_2,\ldots,u_k\).
A survival time of \(t_i=10\) means that subject \(i\) died at time 10. A right-censored time \(u_i=10\) means that at time 10, subject \(i\) was still alive and that we have no further follow-up.
For the moment we will assume that the survival distribution is exponential and that all the subjects have the same parameter \(\lambda\).
We have \(m\) exponential survival times of \(t_1, t_2,\ldots,t_m\) and \(k\) right-censored values at \(u_1,u_2,\ldots,u_k\). The log-likelihood of an observed survival time \(t_i\) is \[
\text{log}\left\{\lambda \text{e}^{-\lambda t_i}\right\} =
\text{log}\left\{\lambda\right\}-\lambda t_i
\] and the likelihood of a censored value is the probability of that outcome (survival greater than \(u_j\)) so the log-likelihood is
Any density on \([0,\infty)\) can be a survival distribution, but the most useful ones are all skew right.
The most frequently used generalization of the exponential is the Weibull.
Other common choices are the gamma, log-normal, log-logistic, Gompertz, inverse Gaussian, and Pareto.
Most of what we do going forward is non-parametric or semi-parametric, but sometimes these parametric distributions provide a useful approach.
1.5 Nonparametric Survival Analysis
Basic ideas
Mostly, we work without a parametric model.
The first task is to estimate a survival function from data listing survival times, and censoring times for censored data.
For example one patient may have relapsed at 10 months. Another might have been followed for 32 months without a relapse having occurred (censored).
The minimum information we need for each patient is a time and a censoring variable which is 1 if the event occurred at the indicated time and 0 if this is a censoring time.
1.6 Example: clinical trial for pediatric acute leukemia
Overview of study
This is from a clinical trial in 1963 for 6-MP treatment vs. placebo for Acute Leukemia in 42 children.
Pairs of children:
matched by remission status at the time of treatment (remstat: 1 = partial, 2 = complete)
randomized to 6-MP (exit times in t2) or placebo (exit times in t1)
Followed until relapse or end of study.
All of the placebo group relapsed, but some of the 6-MP group were censored (which means they were still in remission); indicated by relapse variable (0 = censored, 1 = relapse).
6-MP = 6-Mercaptopurine (Purinethol) is an anti-cancer (“antineoplastic” or “cytotoxic”) chemotherapy drug used currently for Acute lymphoblastic leukemia (ALL). It is classified as an antimetabolite.
Study design
Clinical trial in 1963 for 6-MP treatment vs. placebo for Acute Leukemia in 42 children.
Pairs of children:
matched by remission status at the time of treatment (remstat)
remstat = 1: partial
remstat = 2: complete
randomized to 6-MP (exit time: t2) or placebo (t1).
The estimated hazard in the placebo group is 4.6 times as large as in the 6-MP group (assuming the hazard is constant over time).
1.7 The Kaplan-Meier Product Limit Estimator
Estimating survival in datasets without censoring
Estimating survival in datasets with censoring
For the 6-MP patients, we cannot ignore the censored data because we know that the time to relapse is greater than the censoring time.
For any time \(t\) in months, we know that 6-MP patients with times greater than \(t\) have not relapsed, and those with relapse time less than \(t\) have relapsed, but we don’t know if patients with censored time less than \(t\) have relapsed or not.
The procedure we usually use is the Kaplan-Meier product-limit estimator of the survival function.
The Kaplan-Meier estimator is a step function (like the empirical cdf), which changes value only at the event times, not at the censoring times.
At each event time \(t\), we compute the at-risk group size \(Y\), which is all those observations whose event time or censoring time is at least \(t\).
If \(d\) of the observations have an event time (not a censoring time) of \(t\), then the group of survivors immediately following time \(t\) is reduced by the fraction \[\frac{Y-d}{Y}=1-\frac{d}{Y}\]
Definition 6 (Kaplan-Meier Product-Limit Estimator of Survival Function) If a time-to-event data set contains \(k\) event times \(t_i\), (\(i \in {1:k}\)), where \(n_i\) is the number of individuals at risk at time \(t_i\) and \(d_i\) is the number of events at time \(t_i\), then the Kaplan-Meier Product-Limit Estimator of the survival function is:
Theorem 8 (Kaplan-Meier Estimate with No Censored Observations) If there are no censored data, and there are \(n\) data points, then just after (say) the third event time
where \(\hat F(t)\) is the usual empirical CDF estimate.
Kaplan-Meier curve for drug6mp data
Here is the Kaplan-Meier estimated survival curve for the patients who received 6-MP in the drug6mp dataset (we will see code to produce figures like this one shortly):
Figure 7: Kaplan-Meier Survival Curve for 6-MP Patients
Kaplan-Meier calculations
Let’s compute these estimates and build the chart by hand:
[R code]
library(KMsurv)library(dplyr)data(drug6mp)drug6mp.v2 <- drug6mp |>as_tibble() |>mutate(remstat = remstat |>case_match(1~"partial",2~"complete" ),# renaming to "outcome" while relabeling is just a style choice:outcome = relapse |>case_match(0~"censored",1~"relapsed" ) )km.6mp <- drug6mp.v2 |>summarize(.by = t2,Relapses =sum(outcome =="relapsed"),Censored =sum(outcome =="censored") ) |># here we add a start time row, so the graph starts at time 0:bind_rows(tibble(t2 =0,Relapses =0,Censored =0 ) ) |># sort in time order:arrange(t2) |>mutate(Exiting = Relapses + Censored,`Study Size`=sum(Exiting),Exited =cumsum(Exiting) |> dplyr::lag(default =0),`At Risk`=`Study Size`- Exited,Hazard = Relapses /`At Risk`,`KM Factor`=1- Hazard,`Cumulative Hazard`=cumsum(`Hazard`),`KM Survival Curve`=cumprod(`KM Factor`) )library(pander)pander(km.6mp)
t2
Relapses
Censored
Exiting
Study Size
Exited
At Risk
Hazard
KM Factor
Cumulative Hazard
KM Survival Curve
0
0
0
0
21
0
21
0
1
0
1
6
3
1
4
21
0
21
0.1429
0.8571
0.1429
0.8571
7
1
0
1
21
4
17
0.05882
0.9412
0.2017
0.8067
9
0
1
1
21
5
16
0
1
0.2017
0.8067
10
1
1
2
21
6
15
0.06667
0.9333
0.2683
0.7529
11
0
1
1
21
8
13
0
1
0.2683
0.7529
13
1
0
1
21
9
12
0.08333
0.9167
0.3517
0.6902
16
1
0
1
21
10
11
0.09091
0.9091
0.4426
0.6275
17
0
1
1
21
11
10
0
1
0.4426
0.6275
19
0
1
1
21
12
9
0
1
0.4426
0.6275
20
0
1
1
21
13
8
0
1
0.4426
0.6275
22
1
0
1
21
14
7
0.1429
0.8571
0.5854
0.5378
23
1
0
1
21
15
6
0.1667
0.8333
0.7521
0.4482
25
0
1
1
21
16
5
0
1
0.7521
0.4482
32
0
2
2
21
17
4
0
1
0.7521
0.4482
34
0
1
1
21
19
2
0
1
0.7521
0.4482
35
0
1
1
21
20
1
0
1
0.7521
0.4482
Summary
For the 6-MP patients at time 6 months, there are 21 patients at risk. At \(t=6\) there are 3 relapses and 1 censored observations.
The Kaplan-Meier factor is \((21-3)/21 = 0.857\). The number at risk for the next time (\(t=7\)) is \(21-3-1=17\).
At time 7 months, there are 17 patients at risk. At \(t=7\) there is 1 relapse and 0 censored observations. The Kaplan-Meier factor is \((17-1)/17 = 0.941\). The Kaplan Meier estimate is \(0.857\times0.941=0.807\). The number at risk for the next time (\(t=9\)) is \(17-1=16\).
Now, let’s graph this estimated survival curve using ggplot():
Figure 8: KM curve for 6MP patients, calculated by hand
1.8 Using the survival package in R
We don’t have to do these calculations by hand every time; the survival package and several others have functions available to automate many of these tasks (full list: https://cran.r-project.org/web/views/Survival.html).
The Surv function
To use the survival package, the first step is telling R how to combine the exit time and exit reason (censoring versus event) columns. The Surv() function accomplishes this task.
Now we are going to compare the placebo and 6-MP data. We need to reshape the data to make it usable with the standard survival workflow:
[R code]
library(survival)library(tidyr)drug6mp.v4 <- drug6mp.v3 |>select(pair, remstat, t1, t2, outcome) |># here we are going to change the data from a wide format to long:pivot_longer(cols =c(t1, t2),names_to ="treatment",values_to ="exit_time" ) |>mutate(treatment = treatment |>case_match("t1"~"placebo","t2"~"6-MP" ),outcome =if_else( treatment =="placebo","relapsed", outcome ),surv =Surv(time = exit_time,event = (outcome =="relapsed") ) )
Using this long data format, we can fit a Kaplan-Meier curve for each treatment group simultaneously:
Since \(\text{S}(t)=\text{exp}\left\{-{\Lambda}(t)\right\}\), the Nelson-Aalen cumulative hazard estimate can be converted into an alternate estimate of the survival function:
km_and_na |>ggplot(aes(x = time, y = surv, col = model)) +geom_step() +facet_grid(. ~ strata) +theme_bw() +ylab("S(t) = P(T>=t)") +xlab("Survival time (t, days)") +theme(legend.position ="bottom")
Kaplan-Meier and Nelson-Aalen Survival Function Estimates, stratified by disease group
The Kaplan-Meier and Nelson-Aalen survival estimates are very similar for this dataset.
Copelan, Edward A, James C Biggs, James M Thompson, Pamela Crilley, Jeff Szer, John P Klein, Neena Kapoor, Belinda R Avalos, Isabel Cunningham, and Kerry Atkinson. 1991. “Treatment for Acute Myelocytic Leukemia with Allogeneic Bone Marrow Transplantation Following Preparation with BuCy2.”https://doi.org/10.1182/blood.V78.3.838.838.
Dobson, Annette J, and Adrian G Barnett. 2018. An Introduction to Generalized Linear Models. 4th ed. CRC press. https://doi.org/10.1201/9781315182780.
Kalbfleisch, John D, and Ross L Prentice. 2011. The Statistical Analysis of Failure Time Data. John Wiley & Sons.
Rothman, Kenneth J., Timothy L. Lash, Tyler J. VanderWeele, and Sebastien Haneuse. 2021. Modern Epidemiology. Fourth edition. Philadelphia: Wolters Kluwer.
Vittinghoff, Eric, David V Glidden, Stephen C Shiboski, and Charles E McCulloch. 2012. Regression Methods in Biostatistics: Linear, Logistic, Survival, and Repeated Measures Models. 2nd ed. Springer. https://doi.org/10.1007/978-1-4614-1353-0.